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Abstract 


An algorithm is presented for computing the water-surface profile in steady- state, 
gradually varied flows in tree type irrigation canal networks, given the discharge at the 
root and the control depths at the downstream points. Occurrence of supercritical 
flows in irrigation canal networks is very raLe and hence only the subcritical flow 
conditions are considered The algorithm is based on the principles of (i) decomposing 
the channel network into as small units as possible using the Depth-First -Search 
(DFS) technique, (ii) Solving the smaller units using an appropriate method such as 
fourth-order Runge-Kutta method and (iii) connecting the solutions for the smaller 
units to obtain the final solution for the whole network using the iterative shooting 
method This algorithm js computationally moie efficient by an order of magnitude 
than the direct method using the Newton- Raphson technique. The algorithm does 
not involve the solution of large matrix equations and hence requires less computer 
storage . An example problem for a binary tree type channel network with 42 nodes, 
41 channels and a total of 429 grid points is solved for illustration purpose. The 
complete details of the algorithm are presented in this study 


Chapter 1 


Introduction 


Many hydraulic engineering activities involve the determination of water-surface pro- 
file along the length of the channel for a specified discharge. These computations, 
usually referred to as gradually varied flow (GVF) computations enable safe planning, 
design and operation of open channels. Computation of GVF profiles is required in 
the analysis of problems such as (i) determination of the effect of a hydraulic struc- 
tures on the upstream and downstream channels, (ii) estimation of flood zone and 
(iii) mitigation of losses due to flooding. Solution for steady flow in an open channel 
is also necessary from the point of view of specifying accurate initial conditions for 
unsteady flow computations. The use of unsteady - flow algorithms in a false tran- 
sient approach to determine the steady flow solutions is inefficient computationally. 
Also, the computations may not converge to the proper solution if the transient flow 
algorithm is inconsistent. 

In many situations water-surface profile computations may be required for steady 
flow in open-channel networks. Open-channel networks occur in a braided river in 
mountainous areas and in river system in a delta. Canal networks also form important 
component of an irrigation system. The water from a source gets distributed to field 
through main canal, branch canals, distributaries, minors and water- courses. Many 
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times it is required to determine the water surface profiles in an irrigation canal 
system for a specified discharge entering the system and for the given control depths 
at the down stream ends. 

Dupuit (1848) was probably the first to attempt the computation of GVF profiles . 
Since then many methods have been found for this purpose. They range from graph- 
ical to numerical integration procedures (Chow 1959, Henderson 1966, Subramanya 
1986 and Chaudhry 1993). The numerical integration procedures are currently the 
most popular because (i) they can be adopted to computer application very easily 
and (ii) their range of application is very wide, The earlier numerical methods such 
as the direct step method and the standard step method (Chow 1959 and Henderson 
1966 ) solve the energy equation between two consecutive sections of a channel . The 
later methods such as single-step and predictor-corrector methods (Subramanya 1986 
and Chaudhry 1993) solve the governing ordinary differential equation for GVF using 
an appropriate numerical technique, 

All the numerical integration procedures mentioned in the previous paragraph are 
well suited for single channel applications. Mathematical solution for a canal network 
is not simple and straight forward . The computations usually involve a trial and 
error procedure. For example, consider the problem of flow division at a channel 
junction as shown in Fig. 1.1. The requirement is to determine the water surface 
levels in channels 1, 2 and 3 given the inflowing discharge into the system, i,e. Q\ 
(Q=discharge in m 3 /s) and the controlling depths at the downstream ends, i e, /13 
and /14 (h=flow depth in m). The usual practice of proceeding the computations 
from downstream end to the upstream end for a subcritical flow situation can not 
be adopted because the discharges Q2 and Q$ in channels 2 and 3 are not known 
apriori. Therefore, one of the discharges, say <22 is assumed and the computations 
are carried out for channel 2, there by determining the water surface elevation at 
the junction node 2. Discharge in channel 3 is given by the continuity equation, i e. 

= <2i - (p2- Computations in channel 3 are carried out with this discharge Q3 
and control depth These computations also give a value for the water surface 
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elevation at node 2. Assumed discharge Q 2 is correct if the water surface elevation 
at the junction node 2 determined using computations for channel 2 is same as that 
determined using computations for channel 3. Otherwise, a different value of <52 is 
assumed and the computations are repeated till a correct value of Q 2 is found. The 
above trial and error procedure can be used for a simple channel network such as 
shown in Fig. 1.1 or for a simple island type flow as shown in Fig. 1.2. However, it 
becomes very tedious if it is to be successfully applied to multi island type flows or 
to an irrigation canal network such as shown in Fig. 1.3. 

Significant amount of research has been done regarding the pipe network problem. 
However, not much attention has been paid to the problem of open-channel net- 
works. Considering the problem as analogous to that of the pipe network problem, 
Wylie(1972) developed an algorithm to compute the flow around a group of islands. 
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His approach is based on node continuity and total node energy Like in a pipe 
network problem, this node-and-link method results in a system of non-linear simul- 
taneous equations These equations are solved by the Newton-Raphson technique 
In Wylies’s method, total length of the channel between two nodes is assumed as a 
single reach to calculate loss of energy between two nodes. This assumption is not 
valid who* the lengths of the channels are considerably large. More over, equations 
in Wylie’s work are not in terms of commonly used variables, namely flow depths 
and discharges. Vreugdenhil(1973) developed a similar method for computing the 
flow in a general multi-island system Chaudhry and Shulte(1986) presented similar 
procedure for analyzing steady flows in a channel system having two parallel chan- 
nels Their formulation is in terms of flow depths and discharges, Later, Shulte and 
Chaudhry(1987) extended their method for application to looped channel networks. 
In their method, a channel i m the system is divided into, several reaches Referring 
to the channel reach between sections j and j + 1 of channel 1 (Fig 1.4), the energy 
equation may be written as 


Qh I v — h I 




X t,j 


(S kj + ) (U) 


where 


h— flow depth 
Q = discharge 

g= acceleration due to gravity 
A= cross sectional area 
Z = elevation of channel bottom 
x = distance along the channel length 
S/= friction slope 

Friction slope, 5/ is computed using the Manning’s equation. 

1 AtR* 


( 1 . 2 ) 


where 
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Figure 1 4* Definition sketch for channel reach 
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n=Manning’s n (roughness coefficient) 
R— hydraulic radius of channel section 


The second equation which can be written between the two sections is the continuity 
equation. 

Qi,j+i =: Q\ (1*3) 


Equations (11) and (1.3) can be written for all the reaches of channel i and 
for all the M channels of the system. This results in 2 equations where as 
the number of unknown in the formulation are 2 +■ 1) since there are jV,+i 
nodes in any channel i. Therefore, 2M additional equations aie lequiied foi a unique 
solution. These additional equations are specified by the end conditions and the 
junction conditions. The end conditions are (l) the specific discharge at entrance 
to the system and (ii) the control depths at all down stream points. The junction 
conditions are (i) the continuity equation and (ii) the energy equations. A junction 
of three channels such as shown in Fig. 1.5 results in the following three junction 
equations. 


“ Qi+I.l ~ Q <+2.1 = 0 

(1.4) 

,/v ,+ 1 + Z'.M+i — ^i+i,i + ^i+i.i 

(1.5) 

,N, + i + = ^H-2.1 + Z t +2,l 

(1.6) 


Equations (1.1), (1.3), (1.4), (1.5) and (1 6) written for all nodes, channels and junc- 
tions along with the end conditions constitute the system of equations for computing 
the water surface profile in a channel netwoik Shulte and Chaudhry{1987) solved 
these non-linear simultaneous equations using the Newlon-Raphson technique The 
Newton-Raphson technique requires the computation of Jacobian of the system and 
the determination of its inverse for every iteration step. Shulte and Chaudhry(l987) 
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devised a procedure for numbering the nodes of a parallel channel system, so that the 
resulting nmoufcvtac fs of band width 3M+1 where, M is the channel number of par- 
allel channels. This helps in reducing the computational cost For complex channel 
networks, however, no procedure for arranging of the equations is available to pro- 
duce a Jacobian matrix of minimum width. Even a small irrigation channel network 
shown in Fig. 1 3 results in 340 equations and hence a 340 x 340 Jacobian matrix if 
each channel is divided into nine reaches. The size of the Jacobian increases if the 
number of reaches in each channel increased to obtain better accuracy. Handling such 
large matrices becomes computationally very intensive Shulte and Chaudhry(1987) 
method requires initial estimates for flow depth and discharge at all the nodes. This 
also makes their procedure inefficient for large systems since the convergence will be 
slower if appropriate initial estimates are not chosen. 

In this study, a new method is developed to compute the water surface profile 
for subcritical flow in an irrigation canal network. A binary tiee type of network 
with specified discharge at the upstream end and specified flow depths at all the 
down stream ends is considered The proposed method is based on an iterative 
technique which does not involve the determination of inverse of large matrices. The 
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proposed method is computationally more efficient than the direct method based on 
the Newton-Raphson technique by an order of magnitude and it requires less computer 
storage. 

Formulation of the problem and the algorithm are discussed in chapter-2. Appli- 
cation of the method to an example problem and concluding remarks of the work are 
presented in chapter-3 and chapter-4 respectively 



Chapter 2 

Problem Formulation and 
Solution Algorithm 


Shuite and Chau dhry( 1987) developed a general method for solving open-channel 
network problems which can be used either for a looped network or for a dendritic 
(tree type) network. The objective of this work is to develop a more efficient algorithm 
for solving the network problems by exploiting special properties of a given network. 
Accordingly,rooted binary tree type channel networks without any loops are considered 
in this study. Such networks are common in irrigation systems* For the sake of 
simplicity, it is assumed that (i) discharge from a single source enters the channel 
at upstream end^ (ii) there are no hydraulic structures, with the system (iii) 
only one channel takes off from the main at a section and (iv) the flow is subcriticai 
throughout the system. However, the principles presented in this study can be easily 
extended to other tree type networks with multiple source and intermediate control 
structures. They can also be extended to cases where more than one branch takes 
off from a channel. The occurrence of the supercritical flow in channels is rare and 
requires special treatment. Analysis of such a situation is beyond the scope of the 
present study. 


u 
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Since our proposed method exploits structure of network in solving the system of 
equations, certain concepts fiom Graph Theory have been utilized. Definition of the 
terms used in this work follows. The junction of channels is called a ‘vertex* or 
‘node* The number of channels that are meeting at a node is called the ‘degree* of 
the node The degree of nodes in a binary tree is either three or one. The nodes with 
degree one are called ‘pendant’ nodes except the node where the discharge from the 
source enters the system. Such node is called the ‘root 1 of the system. Referring to 
Fig. 2.1, node 1 is the root of the system and the nodes 4, 5, 8, 9, 11 and 12 are the 
pendant nodes. Each channel is identified by the number of the node downstream of 
the channel. 
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2.1 Basic Governing Equations. 


The objective (s compute water surface elevations (GVF profile) and dischargejthrough- 
out the system given the discharge at the root and the controlling depths at the pen- 
dants. Therefore, the following GVF equation is used to describe the variation of flow 
depth in any channel between the upstream and downstream vertices. 


dh 

dx 


So-Sf 

i _ eay: 

1 5 / 1 3 


(2.1) 


where 


So = slope of the channel bed 
Q = discharge in the channel 
T = water-surface width at any cross section 
/? — momentum coefficient 

The other parameters are as defined earlier (3 is taken as unity in this study since 
we are not concerned with natural irregular channels at this stage . It is to be noted 
here that the standard step method can also be used instead of equation ( 2 . 1 ) for 
computing the water surface profile in a channel between the vertices. Referring to 
the junction formed by the channels i, j and k as shown in Fig 2 2, the continuity 
equation can be written as 

Q, = Qj + Qk ( 2 . 2 ) 

The energy equation can be written as 

Zid + A , d + “—y = Zju + ^ju +■ a 72 " = %ku + hku + 0 -%- (2-3) 

2 gA~ d 2 gA ku 

where the subscript ‘d’ refers to the values at the downstream point of a channel and 
the subscript *u’ refers to the values at the upstream point of the channel. At the 
junction, the form losses and the difference in the velocity heads are generally small 
enough to be considered negligible. Thus neglecting these terms, equation (2.3) may 
be written as 


^id *4* d - fiju Zf; u “h hic U 


(2.4) 
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It is to be noted here that relaxation of the above assumption does not introduce 
complications into the algorithm It can also be seen that any control structure can 
be treated as a junction point with an appropriate governing equation for the energy 
conservation. Equation (2.1) is used for computing the flow depth between two nodes 
while the Eqs. (2.2) and (2 4) are used to determine the flow depth variation at a 
node. 


2.2 Network Algorithm 

The algorithm for successive application of the governing equations to a channel 
network is discussed in this section. The particular network shown in the Fig. 2,1 
is considered for the purpose of explanation As mentioned before, the discharge 
entering the system at the root and the controlling depths at the pendant nodes are 
specified. The channel characteristics are also given for all the channels. It is required 
to compute the water surface levels in all the channels. The first step in the proposed 
method is to assume the discharge at a pendant node, say node 12. With this, the 
discharge and the controlling depth at the downstream end for channel [12] are known. 

i 

\ 


i 
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Therefore, eq. (2.1) can be solved for channel [12] as an initial value problem (IVP) 
by marching the computations from node 12 to the node 10. At the junction node 10 
eq. (2.4) can be applied to determine the upstream depth in channel [11] since the 
upstream depth in channel [12] is known from previous computations. With this, both 
upstream and downstream depths, i.e at nodes 10 and 11 respectively are known 
for channel [11]. However, the discharge is not known. Therefore, eq. (2 1) is solved 
for channel [11] as a Boundary Value Problem (BVP) for a single channel using the 
“shooting method”. This computation determines the discharge in channel [11] as a 
part of its solution. The continuity equation, eq (2.2) and the energy equation, eq. 

(2 4) can now be applied at node 10 for detei mining the discharge and downstream 
depth for channel [10]. Computations can then be marched in channel [10] from the 
downstream node 10 to the upstream node 6 using the IVP subroutine. At node 6, 
eq (2.4) can be applied to determine the upstream depth of channel [7]. However, 
neither the discharge nor the downstream depth m channel [7] are known. Therefore, 
computations at junction node 6 are different from the computations at junction node 
10. The group of channels [7], [8] and [9] are considered as shown in Fig. 2.3. For 
this subsystem, flow depths at the nodes 6,8 and 9 are known* Therefore it can be 
solved as a boundary value problem for a gi'uup of channels. The computations are 
carried out as follows. j 

I 

! 

i 

t 
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• Assume discharge in channel [8] 

• Apply IVP in channel [8] 

• Apply junction conditions at node 7 
t Apply BVP in channel [9] 

• Apply IVP in channel [7] 

• If the computed depth at node 6 is same as the prescribed depth at node 6, the 
assumed discharge in channel [8] is correct Otherwise, it is updated using the 
well known ‘Shooting Method 1 Shooting method is discussed in detail later, 


The above computations give the discharge in channel [7], Equations (2 2) and (2.4) 
can now be applied at node 6 to determine the discharge and downstream depth in 
channel [6], Solution for an IVP can be applied now to march the computations in 
channel [6] from node 6 to node 2 It can be seen from Fig 2 1 that one has to apply 
the nu'lhod (BVP)for group of channels at node 2 in order to determine the discharge 
in channel [3], Discharge in channel [2] is sum of the discharge in channels [6] and 
[3], If this computed discharge in channel [2] is same as the prescribed discharge at 
the root, the assumed value of discharge in channel [12] is correct. Otherwise, it is 
updated using, once again the shooting method. These computations are repeated 
till convergence . The flow chart for the above procedure is shown in Fig. 2 i 


2.2.1 Operation Count 


The efficiency of the proposed algorithm as compaie^to the direct method using 
the N-R technique may be estimated by determining the number of floating point 
operations performed by each method. It is assumed that the forth-order Runge- 
Kutta method used for solving the IVP. The shooting method for the BVP also 
uses the forth-order Runge-Kutta method with iterations. Let each channel in the 
m«-Hoorfe be 'represented 


Junction Eqs at 10 


z 


IVP for [10] to gel h6/10 
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by N grid points, i.e N-l reachs. The determination of water surface profile slope 
at any grid point constitutes one unit operation which equalent to ten floating point 
operations. With this, the application of IVP subroutine to a single channel by the 
fourth order Runge-Kutta method amounts to 40N floating point operations The 
application of BVP subroutine to a single channel amounts to 40pN operations where, 
p is the number of iterations for the convergence of the BVP. Assuming that the same 
number (p) of iterations are taken when the shooting method is applied to a group 
of channels, the total number of floating point operations (refer to the flow chart in 
Fig. 2 4) involved in the proposed method for the channel network of Fig 2 1 , OC p 
is given by 

OC p = 40jV( 1 +3p + 5p 2 + 2p 3 ) (2 5) 

The size of the Jacobian for the channel netwoik shown in Fig. 2 1 is 22N since 
each channel is represented by N grids points and there are two unknowns at each 
grid point. The number of floating point operations involved in the determination of 
the Jacobian is very less compared to the number of operations involved in inverting 
it. The Jacobian is neither banded nor symmetric for a general tree type channel 
network, nor is it diagonally dominant Therefoie, one has to use elimination methods 
for inverting the Jacobian The number of operations for the solution of a 227V x 22A 
matrix equation is equal to (assuming that an LU decomposition method is 

used) Assuming that q iterations are taken for the final convergence of the solution, 
the total operation count for the direct method, OCd js given by 

0C d = 35497V 3 q (2.6) 

The efficiency factor of the proposed method as compaied to the dnect method is 


given by 

OC d 3549 N'q 

(2 7) 


^ ~ 0C P ~ 40N(] +3 p + 5 p 2 + 2P 3 ) 

or 

&9N 2 <i 

~ p ( 3 + 5p + 2 p 2 ) 

(2.8) 


Both the ‘shooting method 1 and the ‘Newton-Raphson’ methods use the gradients 
for determining the update values. Therefore, the number of iterations taken for 
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convergence would be of the same order of magnitude in both the cases Therefore, 
<j> is given by 


0 « 


S9N 2 

3 4 5p 4 2p 2 


(2 9) 


The numerical experience with a targe network in the present study (see chapter 3) 
indicated that p is of the order 10. For a value of p=10, <f> m 0 35 N 2 It can be 
seen that the efficiency factor for the proposed method increases as the number of 
grid points in each channel is increased (j) is equal to 140 foi N =20, which means 
that the proposed method is approximately 140 times more efficient than the direct 
method for the network in Fig. 2 1 , if each channel in the network is divided into 19 
reaches. The proposed method does not involve handling of large matrices and hence 
the memory requirements are also small 


2.3 Path of Marching 


The proposed algorithm results in high efficiency when the calculations march along 
a certain path in the network For example, consider the case of staiting the compu- 
tations at the pendant node 5 of the example netwoik (Fig 2 1) The flow chart of 
the application of the proposed algoiithm is shown in Fig. 2 5. The operation count 
for this case is equal to 40jV(1 4 2p 4 4p 2 4 3p 3 4 p 4 ) The efficiency factor, <f> is 
approximately equal to 0 066 Af 2 , which is 5 times less when compared to the applica- 
tion of the proposed algorithm along the path given in Fig. 2 4 . Theiefoie, the “best 
path of marching” has to be determined before starting the computations, The best 
path of marching determines (i) the pendant node at which the computations begin 
for the whole network and (ii) the pendant node at which the computations begin 
while solving the BVP for a group of channels. 
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Figure 2 6 Definition sketch foi node classification 

2.3.1 Node Classification 

It is required to represent the network as a ‘graph’ in order to determine the best 
path of marching. In a binary tree type network such as considered in the present 
study, a node receives water from only one node and distiibutes the water to two 
other nodes (Fig 2.6) . The node from wheie watei is received, J v is termed as the 
‘parent’ of node J and the nodes to where water is distributed, J\ c and J rc aie termed 
as ‘children’ of the node J, One of them, Ji c is teimed as the ‘left child’ and the 
other J re is termed as the ‘right child’ We can easily represent the graph in terms of 
numbers specifying the parent, left child and right child of all the nodes in the graph 
Zero is specified if a node does not have either a paient, a left child or a light child. 
A node with zero parent is the ‘root’ A node with no children is the ‘pendant 5 . A 
node which has a pendant node either as a left child or as a right child or as both is 
termed as Type- 1 node. A node which does not have a pendant node as a left or as a 
right child is termed as Type-II node. Referring to Fig. 2.7. nodes 5, 7, 8, 10, 11, 13, 
17, 18, 19, 22, 23, 25, 27 and 28 are pendant nodes; nodes 4, 6, 9, 12, 15, 16, 21, 24 
and 26 are Type-I nodes and nodes 2, 3, 14 and 20 are TypeTI nodes. It can be seen 
that the the computations at a type-I node involve the solution of BVP foi a single 
channel, where as the computations at a Type-II node involve the solution of BVP 


CENTRAL LIBRARY 

' I l. T„ KANPUa 




23 


for a group of channels 


2.3.2 Algorithm for Path of Marching 


Consider the computations at node 2 of the network. These computations involve 
matching the depth, h 2 / 3 (depth at node 2 in channel [3]) with the depth, k 2 /i 2 ~ 
(depth at node 2 in channel [12]) This can be done m two ways In alternative 1, 
h 2 /3 is determined first and then the BVP foi gioup of channels on the light side 
of 2 is solved to determine /12/12 In the alternative 2, h 2 /] 2 is detei mined hist and 
then the BVP for group of channels on left side of 2 is solved to determine A 2 /3 The 
BVP for left side of 2 involve the BVP for group of channels at 3 and the BVP for 
single channel at 4, 6 and 9, i.e one BVPGC (BVP for group channel) and three 
BVPSC’s(BVP for single channel). On the other hand , the BVP for the right side 
of 2 involves two BVPGC’s at nodes 14 and 20 and six BVPSC’s at nodes 12, 15, 

16, 21, 24 and 26 . This means, the BVP for left side of 2 involves less computations 
than the BVP for right side of 2. Therefore it is better to make computations for the 
right side of 2 first^then iterate on the left side, we alternative 2 should be chosen. 
Similarly, while doing computations for node 14, it is better to do the computations 
on the left side first and then iterate for BVPGC on the right side since there is no 
type-II node on the right side. There are no type- II nodes downstream of node 20 
either on the left side or on the right side. However, the number of type-I nodes 
(two ) on the left side is greater than the numbei of type-I nodes (one) on the right 
side. Therefore, while doing computations for node 20, it is bettei to peiform the 
computations on the left side first and then iterate for BVPGC on the right side 
Similarly for node 3, the computations on the light side aie peifoimed fiist and then 
the iterations are performed for the BVPGC on the left side. With the above, the 
path of marching as shown in Fig 2 S is obtained Best path of matching is found 
such that a minimum number of BVP’s (see Fig 2 8) occui inside the nested loops. 
The sequence for computations is based on the fact that computations at a type- 
II node are more intensive than the computations at a type-I node. The following 
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Figure 2.8: Flow chart for network in Fig. 2.7 
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rules are established for finding the best path of marching for computation. Channel 
network shown in Fig. 2.7 is considered for lllustiation purposes 

1 At every node, the number of type-II nodes downstream of the node on the 
left side, N/f/i and the number of type-II nodes on the right side, ./V////? 
are counted The number of type -I nodes on the left side N/jl and the 
number of type- 1 nodes on the right side Nijr are also counted 

2. At a type-II node, if Nu/n > N }J / L then the left side is the iterating side 
and the right side is the direct computing side and vice versa For example, 
right side is the iterating side at node 14 

3 At a type-II node, if Nii/r = Nii/i then the number type-I nodes are 
compared If Ni/r > Ab/z, then the direct computations are made for the 
right side and the iterations are performed on the left side and vice versa. 
For example, right side is the iterative side at node 20. 

4. AT a type-I node, the branch connecting the pendant node is always solved 
as a BVP The direct computations are made for the other side For 
example, BVP is applied for channel [13] 

5. The objective is to choose the flow in a pendant channel for performing 
iterations at a type-II node such that the same flow variable is not chosen 
as the iterative variable for iterations at two different type-II nodes This 
is achieved by adopting the following procedure. 

The downstream channels at every node are labeled either as iterative 
channel or as direct computation channel. There are two pendant nodes 
reachable from every type-II node along selected paths, one on the direct 
computation side and the other on the iterative side, denoted as D- node 
and I-node respectively. The D-node coi responding to a type-II node is 
reached by following the direct computation channels at every node along 
the path. The I-node corresponding to a type-II node is leached by first 
moving along the iterative channel of the type-II node and subsequently 
following the direct computation channel at every other node along the 
path. 
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6. The computations for the whole network always stall at the D-node cor- 
responding to the first node closest to the loot Discharge in this D-node 
pendant channel is assumed and iterated upon as the outer loop. For ex- 
ample, computations for the channel network in Fig 2 7 start with the 
assumption for discharge in channel [27). 

7. The iterations at a type-II node start at the corresponding I-node. For ex- 
ample iterations at node 14 start with an assumption for Q i 8 and iterations 
at 2 start with an assumption for Q 8 . 

8. The BVP for single channel is applied only foi those pendant channels 
which are not chosen as iterative vanables and is never applied for any 
intermediate channel. The IVP for an intei mediate channel can be applied 
only after flow depths in both the downstream channels at the junction 
are matched. For example, IVP for channel [20j can be applied only after 
obtaining ^ 20/24 = ^ 20 / 21 * 


The above principles are used to develop a progiam for determining the path of com- 
putation and the sequence for calling the IVP and BVP subroutines for a given channel 
network. Depth first search (DFS) technique is used foi developing the piogram (Deo 
1974). 


2.4 Initial Value Problem 


In the context of GVF computations foi a single open-channel flow, the initial value 
problem(IVP) refers to the solution of eq. (2.1) for a specified discharge, Q and 
the given downstream flow depth, h d (Fig 2 9). This can be easily solved using 
fourth-order Runge-Kutta method (Subramanya 1987, Chaudhry 1993) 
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h ^{specified) 


Channel bed 



2.5 Boundary Value Problem 


The BVP for a single channel refers to the solution of differential equation, eq (2.1) 
for specified upstream and downstieam flow depths, h u and hj respectively (Fig. 
2.10) The BVP is solved using the “Shooting Method". In this method, the BVP is 
solved as an IVP with iterations The discharge in the channel, Q i is assumed first 
The IVP for the channel is solved using the specified downstieam depth and the 
estimated discharge to obtain the upstream flow depth, A u ,j. h Ut i is compared with 
the specified flow depth, A u If k Ut i is less than h u then the discharge in the channel 
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is increased by a small amount^ Glse, the dischaige is decreased. A change of b% 
is used in the present study. /i Uj2 , conesponding to the dischaige, Q 2 is delei mined 
using the I VP subroutine again The next approximation Qz is determined by a 
simple extrapolation as given below. 

<?3 = Qi + (Q, - Qj) ( 2 . 10 ) 

"u,l — n u ,2 

Q 2 is then assigned to^Qi and Q 3 is assigned to Q 2 . The new value of h u2 corre- 
sponding to the new value of Q 2 is determined using the IVP and the procedure is 
repeated till convergence. The convergence criteria is, 

\h u - /i Ui2 | < e (2.11) 


in which t is the specified tolerance 


2.6 Iteration at a Type-II Node 


It is required to solve the BVP for a group of channels at every type-II node. For 
example, consider the iterations at node 14 of the channel network shown in Fig 2.7 
, The corresponding I-node for node 14 is node 18 The iterations at node 14 start 
with an assumption for the discharge in channel [18]. The flow depth at node 14 in 
channel [15], i.e. / 114/15 is determined with the following sequence of operations. 


• Call IVP for [ 18 ) 

• Call BVP for [17] 

• Apply continuity eq. at node 16 

• Call IVP for [16] 

• Call BVP for [19] ' 
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• Apply continuity eq. pinocle 

• Call IVP for (15) 


The above determined value of ^ rh* compared with h\^j 2 o (flow depth at node 
14 in channel [20]}, cohieh lias l>ee^ wi^jputed at a previous level - Discharge 
is revised if kufu is not equal Uo „ The principles of “Shooting Method” as 

explained in the previous section for this purpose The same procedure is 

followed at all the type-II nodes inajy network (Refer to Fig 2.6). 


2.7 Iterations for the (Puter Loop 


The iterations for the outer loop to the iterations at the first node which is 
closest to the root. The i teralious^Wi the corresponding D-node with an initial 
estimate for the discharge in that j^inLaAt channel The computations are performed 
as explained before to deitermin^ th^iigtharge m the channel connected to the root 
Estimated discharge in iteration, vhatfbel is revised if the above computed discharge 

is not equal to the specified dis«chu,^ gt the root. Shooting method is used for this 
purpose, 


2.8 Initial Estimate 


Initial estimates of discharges r^T) il'-ed for staiting the calculations as well as for 
obtaining the solution of GvPdETtyeS and Type-II nodes. This means the initial 
estimates of discharges are iw] uii^d forf all the pendant channels. In this study, the 
initial estimate discharge in any p^nd^bt channel is approximately equal to the total 
discharge divided by the nutnkerAf pedant nodes. However, based on the previous 
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computations, the initial estimates of dischaiges foi obtaining the solution of BVP 
can be revised as the outer loop iterations proceed 


2.9 Conditioning of the Iterations 

The proposed algorithm is based on two assumptions, (i) the flow is subcntical 
throughout the network and (ii) the flow direction is always from the loot towards 
the pendants. Both these assumptions are valid for the noimal steady slate operating 
conditions in irrigation canal network. For a given network and the end conditions, 
the actual solution may satisfy these two conditions But, the flow may either become 
negative or may become supercritical during intermediate iterations while applying 
the proposed method as outlined in the previous sections Thi« happens when the 
initial estimates for the discharges are for from correct solution Therefore, internal 
numerical controls are built into the proceduie which enable the computations to 
march forward even when the above situations arise. 


2.9.1 Occurrence of Negative Flow 


Consider the computations at a type-II node, foi example take node 14 of Fig 2.7, 
corresponding to an outer loop iteration. The calculations up to this stage of outer 
loop iterations have determined the depth at this type-II node on thediiect computing 
side, i.e. hu /20 Inner loop iterations have to be performed at this node to mahct with 
the depth on the direct computing side, i.e h u / 2 o + Zh /2 0 = < u /15 + i/is This 
amounts to solving the boundary value problem foi gioup of channels with depths 
at the type-II and all the pendant nodes, i.e nodes 17, IS and 19 being specified 
The computation start with an initial estimate Toi the discharge in the <.lia.ni.jcl 
connecting the corresponding I-node, i e. initial estimate for Q is The iterations using 
shooting method are performed to update this value of Q 18 . However, the discharge 


31 


Qia rnay come out to be negative during these iterations. This could be because the 
depth at type-II node, i . e. as determined in the previous calculations is very 

small compared to the control depth at the corresponding 1-node, i.e. /i^g. When 
this happens, the discharge in iteration channel, i.e Q \ g is set to zero, and one more 
iteration of the BVP is performed. If the the updated discharge remains negative, then 
the iterations for the BVP are terminated even though the convergence criterion is not 
satisfied. The continuity equation is applied at the type-II node, i e. node 14 using 
the most curient values available for the flow discharges and the computations are 
marched further. A similar technique is adopted when the discharge becomes negative 
during iterations for BVP at type-I node. The above procedure for conditioning the 
iterations is based on the logic that (i) the flow in any channel is always from parent 
to a. child node and (ii) the convergence of the iterations for an inner loop such as 
the iterations of BVP at a type-II or type-I node is strictly required only for the 
final outer loop iteiation. Although this conditioning results m spurious discharge 
values for intermediate iterations, this situation does not occur as the solution from 
the outer loop converges to the correct value 


2.9.2 Occurrence of Supercritical Flow 


The flow may become supercritical either at a node or in a channel during an in- 
termediate iteration. This situation may arise either while performing an outer loop 
iteration or while performing iterations at a type-II or type-I node. Considering the 
iterations at a type-II node, supercritical flow condition may occur if the upstream 
flow depth for the BVP is very large compared to the specified depths at the cor- 
responding pendant nodes Similar situation may arise while solving the BVP for a 
single channel if the upstream and downstream flow depths are incompatible from the 
point of view of sustaining the subcritical flow throughout Supercritical flow condi- 
tions may arise while solving the IVP for an intermediate channel if the water surface 
profile for that discharge and downstream depth happens to be of the type SI. In the 
proposed method , the Froude number of the flow is calculated at every grid point of 
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the channel as pari of computations for the water surface profile. When the Froude 
number becomes greater than 0.95, the discharge in the channel is reduced in a step- 
wise manner till the flow becomes subcritical for the entire channel This procedure 
may not satisfy the convergence criterion at the upstream node if the supercritical 
flow occurs while solving the BVP for a single channel. In such a case the discharge 
is updated as before and iterations are continued If the flow becomes supercritical 
in two consecutive iterations, then the iterations for -this type-I node are terminated 
even though the convergence criterion is not satisfied The continuity equation is ap- 
plied using the most current values available for the discharge and the computations 
are marched ahead. The computations are marched beyond the upstream node in 
the normal way if the problem of supercritical flow occurs while solving the IVP for a 
channel. The above procedure for conditioning the iteration affect the convergence of 
iterations at a Type- II node if the TVP or the BVP for a single channel fall within the 
iterative loop for this node The only possible remedy for this problem is to apply a 
constraint for maximum number of iterations The above procedure for conditioning 
the iterations is based on a similar logic as explained earlier, i e. the convergence of 
iterations for an inner loop is strictly required only for the final outer loop calculation 


2.10 Closure 

A computer program incorporating the principles presented in the previous sections 
has been developed for computing the steady state GVF profiles in binary tree type 
irrigation canal networks The inputs for executing this program are (i) Channel char- 
acteristics and (ii) the node connectivities. The program first sorts out the network 
for the node classification and determines the path of computations for the outer loop 
as well as for the iterations at the interior points. According to the above path of 
computation, the program then accesses the subroutines for the initial value problem 
and the boundary value problem in a sequential way to determine the water surface 
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profile The application of the algorithm to an example network is presented in the 
next chapter. 



Chapter 3 


Illustrative Example 


In this chapter, an idealized channel network shown in Fig 3.1 is analyzed with the 
computational procedure described m the previous chapter. The nodes are numbered 
as shown in that figure and the arrows indicate the flow directions. All channels 
have trapezoidal cross-sections. Form losses and the difference in velocity heads at 
junctions are assumed to be zero. The bed levels on all the three sides of a junction 
are also assumed to be same. The flow in all the channels is subcritical , and the 
end conditions are the specified discharge in channel [2], Qi = 40 m 3 /s and the 
flow depths at all the pendant nodes as given in table 3.1 The data for channel 
characteristics are listed in table 3.2. 

The computations start with identifying the nodes into type-I and type-II. The 
node type is indicated for all the nodes in Fig 3d. The I-nodes and D- nodes corre- 
sponding to all the type-I I nodes obtained using the proposed method are shown m 
table 3 3. The sequence of accessing the IVP and BVP subroutines and the loops for 
the iterations are shown in Fig. 3 2. The iterations for surface profile computation 
start with the assumption for discharge in channel [28] since node 28 is the D-node 
corresponding to node 2. The initial estimate for Q 28 =2.0 m 3 /s, which is nearly 
equal to the total discharge divided by the total number of pendant channels. The 
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Figure 3.2: Flow chart for example network 
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initial estimates for discharges at all the I-nodes, i.e. in all iterating pendant channels 
as well as the initial estimates for the discharges in all the BVPs for single channels 
are also equal to 2.0 m 3 /s. A tolerance of 5 x 10“ 4 m is used as the convergence 
criteria for the BVP for single channels as well as the convergence of iterations at a 
type-II node. A tolerance of 0-01 m 3 /s is used for the convergence of iterations for the 
outerloop. This corresponds to an average error of 5 x 10“ 4 m 3 /s for the discharge in 
a pendant channel. The discharge, the upstream flow depth and the downstream flow 
depth for all the channels obtained using the proposed method are shown in table 
3.4. The values given in table 3*4 are verified by computing the flow depths at all the 
nodes except the root using the forward problem. The discharges in all the channels 
and the flow depth at the root already calculated by the proposed method are used 
directly in the forward problem. The flow depths computed by the forward problem 
are identical (within the specified tolerance) with the values presented in table 3.4. 

The method converged to the correct solution within six iterations of the outer loop. 
The number of iterations performed at any type-II node varied from a maximum value 
of 9 to a minimum value of 3. The number of iterations performed for any BVP for 
a single channel varied form a maximum value of 7 to minimum of 2. The number of 
iterations at any type-II or type-I node decreased as the solution from the outerloop 
converged to the correct value. Assuming that the proposed method takes seven 
iterations for the outerloop, nine iterations at every type-II node and seven iterations 
for the BVP for a single channel, the operation count for the proposed method is 
approximately equal to 6.8 X 10 6 . The operation count for the application of the 
direct method with seven iterations would be equal to 1.5 x 1G 9 . Therefore, the 
proposed method is approximately 220 times more efficient than the direct method 
for solving the example network problem. The actual efficiency is higher than the 
above since the number of iterations at any type-I or Type-II node decreases as the 
solution for the outerloop converges to the correct value. The other advantages of the 
proposed method are, (i) It does not involve the solution of large matrices and hence 
the memory requirements are very low and (ii) The initial estimates are required for 
the discharges in the pendant channels only. 


i 
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Table 3.1: Flow depths at pendant, nodes 


Node 

Depth(m) 

Node 

Depth(m) 

Node 

Depth(m) 

5 

0.9111 

25 

1.3622 

36 

1 7769 

9 

1.6559 

28 

1.4766 

37 

1.2190 

12 

0.9759 

30 

1.1741 

38 

1.4745 

15 

0.9127 

32 

1.0749 

39 

1.3719 

18 

1.6021 

33 

1,4777 

40 

1.6091 

20 

1.8784 

34 

1.7107 

41 

13310 

22 

1.6729 

35 

2.0070 

42 

1.2535 
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Table 3 . 2 : Channel characteristics for network in Fig . 3.1 


Channel Lenth ( m ) Bed width ( m ) Side slope Bed sloped n 




10.0 



8.5 

4 

1700.0 

7.0 

5 

1500.0 

5.0 

6 

1500 0 

5.0 

7 

1400.0 

4.0 

8 

1200.0 

3.0 

9 

1000.0 

2.0 

10 

1400.0 

3.5 

11 

1200.0 

2.7 

12 

1000.0 

1.75 

13 

1300.0 

25 

14 

1200.0 

1.5 

15 

1000.0 

1.0 

16 

1000.0 

1.5 

17 

1000.0 

1.0 

18 

1000.0 

1.75 

19 

1000.0 

1.5 

20 

900.0 

0.9 

21 

1100.0 

1,5 

22 

1000.0 

1.0 

23 

1200.0 

1.75 

24 

1100.0 

1.5 

25 

1000.0 

1.0 

26 

1200.0 

2.0 

27 

1000.0 

1.75 

28 

900.0 

1.5 

29 

900.0 

1.5 

30 

800.0 

1.0 

31 

800.0 

1.25 

32 

700.0 

0,75 

33 

700.0 

0.5 

34 

700.0 

0.5 

35 

700.0 

0.5 

36 

700.0 

0.5 

37 

700.0 

0.5 

38 

700.0 

0.5 

39 

700.0 

0.5 

40 

700.0 

0.5 

41 

700.0 

0,5 

42 

700.0 

0.5 


2.0 

7500.0 


2 0 

6667.0 


2.0 

6250.0 

0.017 

2.0 

5882.0 

0.018 

2.0 

5000.0 

0.020 

2.0 

4762.0 

0.020 

2.0 

4545.0 

0.020 

2.0 

4166.0 

0 022 

1.0 

4000 0 

0.022 

1.0 

4545 0 

0.022 

2.0 

4166 0 

0.022 

20 

4545.0 

0.022 

1.0 

4000 0 

0.022 

20 

4545.0 

0.022 

2.0 

4166.0 

0.022 

1.0 

4000.0 

0.022 

2.0 

4166 0 

0 022 

20 

4166.0 

0.022 

0.9 

4000 0 

0.022 

2.0 

4166.0 

0.022 

1.0 

4000.0 

0.022 

2.0 

4166 0 

0.022 

2.0 

4166.0 

0.022 

1.0 

4000.0 

0 025 

20 

4166.0 

0.02 

2.0 

4166.0 

0.022 

2.0 

4166,0 

0.022 

1.0 

4000.0 

0.022 

1.0 

4000.0 

0.022 

2.0 

4166 0 

0 022 

2.0 

4166.0 

0.022 

1.0 

2000.0 

0.03 

1.0 

2000.0 

0 030 

1.0 

2000.0 

0.03 

1,0 

2000.0 

0.03 

1.0 

2000.0 

0.03 

1.0 

2000.0 

0.03 

1.0 

2000.0 

0.03 

1.0 

2000.0 

0.03 

1.0 

2000.0 

0.03 

1.0 

2000.0 

0,03 


No.of reaches ( N ) 
20 
20 
15 
15 
15 
15 
15 
10 
15 
15 
15 
15 
15 
15 
10 
10 
10 
10 
10 
10 
10 
10 
10 
10 
10 
10 
10 
10 
8 
8 
8 
5 
5 
5 
5 
5 
5 
5 
5 
5 
5 


I 


» 40 


I 

i 
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Table 3.3: I-nodes and D-nodes for type-II nodes 


Type-II node 

I-node 

D-node 

2 

35 

28 

3 

32 

28 

6 

18 

35 

7 

22 

35 

10 

38 

28 

13 

41 

32 




Table 3.4: Solution for example network in Fig. 3.1 


Channel 

Discharge 

(Q) 

Upstream 
depth (in) 

Downstream 
depth (tyO 

2 

39.9958 

2.1921 

1.8664 

3 

24.6298 

1.8664 

1.6821 

4 

16.9722 

1.6821 

1.5045 

5 

11.3176 

1.5045 

.911099 

6 

15.366 

1.8664 

1 6503 

7 

10.0897 

1.6503 

1.5165 

8 

6.352 

1.5165 

1.5333 

9 

3.8329 

1.5333 

1 6559 

10 

7.6576 

1.6821 

1 3403 

11 

3 9379 

1.3403 

1 0036 

12 

1.9537 

1.0036 

975899 

13 

5 6547 

1.5045 

1.4048 

14 

2 8783 

1.4048 

.9845 

15 

1.4032 

.9845 

912699 

16 

5.2762 

1.6503 

1.6352 

17 

3.0786 

1.6352 

1.4232 

18 

2.3152 

1.4232 

1.6021 

19 

2,5191 

1.5333 

1 7128 

20 

1.4325 

1.7128 

1.8784 

21 

3.7377 

1.5165 

1.6025 

22 

2.2065 

1.6025 

1.6729 

23 

3.7197 

1.3403 

1 2658 

24 

2.5728 

1 2658 

1.3424 

25 

1.4612 

1.3424 

1.3622 

26 

1.9842 

1.0036 

1.1098 

27 

1.4624 

1.1098 

1.2815 

28 

1.0853 

1.2815 

1.4766 

29 

2.7764 

1.4048 

1.2693 

30 

1.7083 

1.2693 

1.1741 

31 

1.4751 

.9845 

1.0016 

32 

1.0338 

1.0016 

1.0749 

33 

2.1976 

1.6352 

1.4777 

34 

.763411 

1.4233 

1 7107 

35 

1.0866 

1.7128 

2.007 

36 

1.5313 

1.6025 

1.7769 

37 

1.1469 

1.2658 

1.219 

38 

1.1116 

1.3424 

1.4745 

39 

.521841 

1.1098 

1.3719 

40 

.377035 

1.2815 

1.6091 

41 

1.0681 

1,2693 

1.331 

' 42 

.441318 

1.0016 

1.2535 



Chapter 4 


Concluding Remarks 


In this study, an iterative method has been developed for computing water surface 
profiles in tree type channel networks. The method solves for the water surface profile 
in a network given the inflow discharge at the root and the control depths at all the 
clown stream points. The proposed method has been used for solving the gradually 
varied flow in a binary tree type channel network with inflow from a single source. 
However, the method can be easily extended to other tree type channel networks and 
to cases where inflow into the network is at several locations. The proposed method 
is very useful for analyzing steady gradually varied flow in irrigation canal networks. 

The proposed method is based on the principles of (i) decomposing the channel 
network into as small units as possible, (ii) solving the small units using an appropriate 
algorithm and (iii) connecting the solution for the small units to obtain the final 
solution using an iterative technique. Each channel in the network is either solved 
as an initial value problem (IVP) or a boundary value problem (BVP) Depth first 
search (DPS) algorithm is first used to determine whether a channel belongs to the 
IVP category or the BVP category. Accordingly, the IVP and BVP subroutines are 
called sequentially to determine the water surface profile. Fourth-order Runge-KuUa 
method is used for solving the IVP. Shooting method is used for solving the BVP as 
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well as revising the intermediate solution during iterations 

For tree type channel networks, the proposed method is more efficient than the 
direct method based on the Newton-Raphson technique (Shulte and Chaudhry, 1987) 
by an order of magnitude. It does not involve solving large matrix equations and 
hence requires less computer storage. Suitable internal controls can be built into 
algorithm to make it. more robust and efficient. However, at the present stage of 
development, the proposed method can not be used for analyzing general looped 
channel networks for which the direct method developed by Shulte and Chaudhry 
(1987) can be applied. 
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